      REAL*8 FUNCTION PRAXIS(T0,MACHEP,H0,N,PRIN,X,F,FMIN)
C                             LAST MODIFIED 3/1/73
      REAL*8 T0,MACHEP,H0,X(N),F,FMIN
      EXTERNAL F
      INTEGER PRIN
C     RETRIEVED FROM:
C     http://www.netlib.org/opt/praxis
C     2 SEP 2009
C
C     PRAXIS RETURNS THE MINIMUM OF THE FUNCTION F(X,N) OF N VARIABLES
C     USING THE PRINCIPAL AXIS METHOD.  THE GRADIENT OF THE FUNCTION IS
C     NOT REQUIRED.
C
C     FOR A DESCRIPTION OF THE ALGORITHM, SEE CHAPTER SEVEN OF
C     "ALGORITHMS FOR FINDING ZEROS AND EXTREMA OF FUNCTIONS WITHOUT
C     CALCULATING DERIVATIVES" BY RICHARD P BRENT.
C
C     THE PARAMETERS ARE:
C     T0       IS A TOLERANCE.  PRAXIS ATTEMPTS TO RETURN PRAXIS=F(X)
C              SUCH THAT IF X0 IS THE TRUE LOCAL MINIMUM NEAR X, THEN
C              NORM(X-X0) < T0 + SQUAREROOT(MACHEP)*NORM(X).
C     MACHEP   IS THE MACHINE PRECISION, THE SMALLEST NUMBER SUCH THAT
C              1 + MACHEP > 1.  MACHEP SHOULD BE 16.**-13 (ABOUT
C              2.22D-16) FOR REAL*8 ARITHMETIC ON THE IBM 360.
C     H0       IS THE MAXIMUM STEP SIZE.  H0 SHOULD BE SET TO ABOUT THE
C              MAXIMUM DISTANCE FROM THE INITIAL GUESS TO THE MINIMUM.
C              (IF H0 IS SET TOO LARGE OR TOO SMALL, THE INITIAL RATE OF
C              CONVERGENCE MAY BE SLOW.)
C     N        (AT LEAST TWO) IS THE NUMBER OF VARIABLES UPON WHICH
C              THE FUNCTION DEPENDS.
C     PRIN     CONTROLS THE PRINTING OF INTERMEDIATE RESULTS.
C              IF PRIN=0, NOTHING IS PRINTED.
C              IF PRIN=1, F IS PRINTED AFTER EVERY N+1 OR N+2 LINEAR
C              MINIMIZATIONS.  FINAL X IS PRINTED, BUT INTERMEDIATE X IS
C              PRINTED ONLY IF N IS AT MOST 4.
C              IF PRIN=2, THE SCALE FACTORS AND THE PRINCIPAL VALUES OF
C              THE APPROXIMATING QUADRATIC FORM ARE ALSO PRINTED.
C              IF PRIN=3, X IS ALSO PRINTED AFTER EVERY FEW LINEAR
C              MINIMIZATIONS.
C              IF PRIN=4, THE PRINCIPAL VECTORS OF THE APPROXIMATING
C              QUADRATIC FORM ARE ALSO PRINTED.
C     X        IS AN ARRAY CONTAINING ON ENTRY A GUESS OF THE POINT OF
C              MINIMUM, ON RETURN THE ESTIMATED POINT OF MINIMUM.
C     F(X,N)   IS THE FUNCTION TO BE MINIMIZED.  F SHOULD BE A REAL*8
C              FUNCTION DECLARED EXTERNAL IN THE CALLING PROGRAM.
C     FMIN     IS AN ESTIMATE OF THE MINIMUM, USED ONLY IN PRINTING
C              INTERMEDIATE RESULTS.
C     THE APPROXIMATING QUADRATIC FORM IS
C              Q(X') = F(X,N) + (1/2) * (X'-X)-TRANSPOSE * A * (X'-X)
C     WHERE X IS THE BEST ESTIMATE OF THE MINIMUM AND A IS
C              INVERSE(V-TRANSPOSE) * D * INVERSE(V)
C     (V(*,*) IS THE MATRIX OF SEARCH DIRECTIONS; D(*) IS THE ARRAY
C     OF SECOND DIFFERENCES).  IF F HAS CONTINUOUS SECOND DERIVATIVES
C     NEAR X0, A WILL TEND TO THE HESSIAN OF F AT X0 AS X APPROACHES X0.
C
C     IT IS ASSUMED THAT ON FLOATING-POINT UNDERFLOW THE RESULT IS SET
C     TO ZERO.
C     THE USER SHOULD OBSERVE THE COMMENT ON HEURISTIC NUMBERS AFTER
C     THE INITIALIZATION OF MACHINE DEPENDENT NUMBERS.
C
      LOGICAL ILLC
      INTEGER NL,NF,KL,KT,KTM
      REAL*8 S,SL,DN,DMIN,FX,F1,LDS,LDT,T,H,SF,DF,QF1,QD0,QD1,QA,QB,QC
      REAL*8 M2,M4,SMALL,VSMALL,LARGE,VLARGE,SCBD,LDFAC,T2,DNI,VALUE
      REAL*8 RANDOM,DSQRT,DABS,RANF
C
C.....IF N>20 OR IF N<20 AND YOU NEED MORE SPACE, CHANGE '20' TO THE
C     LARGEST VALUE OF N IN THE NEXT CARD, IN THE CARD 'IDIM=20', AND
C     IN THE DIMENSION STATEMENTS IN SUBROUTINES MINFIT,MIN,FLIN,QUAD.
C     MODIFIED TO ACCOMODATE 100 FREE PARAMETERS
C     2 SEP 2009
C
      REAL*8 D(100),Y(100),Z(100),Q0(100),Q1(100),V(100,100)
      COMMON /GLOBAL/ FX,LDT,DMIN,NF,NL
     .       /Q/ V,Q0,Q1,QA,QB,QC,QD0,QD1,QF1
C
C.....INITIALIZATION.....
C     MACHINE DEPENDENT NUMBERS:
C
      SMALL=MACHEP*MACHEP
      VSMALL=SMALL*SMALL
      LARGE=1.D0/SMALL
      VLARGE=1.D0/VSMALL
      M2=DSQRT(MACHEP)
      M4=DSQRT(M2)
C
C     HEURISTIC NUMBERS:
C     IF THE AXES MAY BE BADLY SCALED (WHICH IS TO BE AVOIDED IF
C     POSSIBLE), THEN SET SCBD=10.  OTHERWISE SET SCBD=1.
C     IF THE PROBLEM IS KNOWN TO BE ILL-CONDITIONED, SET ILLC=TRUE.
C     OTHERWISE SET ILLC=FALSE.
C     KTM IS THE NUMBER OF ITERATIONS WITHOUT IMPROVEMENT BEFORE THE
C     ALGORITHM TERMINATES.  KTM=4 IS VERY CAUTIOUS; USUALLY KTM=1
C     IS SATISFACTORY.
C
      SCBD=10.D0
      ILLC=.FALSE.
      KTM=1
C
      LDFAC=0.01D0
      IF (ILLC) LDFAC=0.1D0
      KT=0
      NL=0
      NF=1
      FX=F(X,N)
      QF1=FX
      T=SMALL+DABS(T0)
      T2=T
      DMIN=SMALL
      H=H0
      IF (H.LT.100*T) H=100*T
      LDT=H
C.....THE FIRST SET OF SEARCH DIRECTIONS V IS THE IDENTITY MATRIX.....
      DO 20 I=1,N
           DO 10 J=1,N
10              V(I,J)=0.0D0
20         V(I,I)=1.D0
      D(1)=0.D0
      QD0=0.D0
      DO 30 I=1,N
           Q0(I)=X(I)
30         Q1(I)=X(I)
      IF (PRIN.GT.0) CALL PRINT(N,X,PRIN,FMIN)
C
C.....THE MAIN LOOP STARTS HERE.....
40    SF=D(1)
      D(1)=0.D0
      S=0.D0
C
C.....MINIMIZE ALONG THE FIRST DIRECTION V(*,1).
C     FX MUST BE PASSED TO MIN BY VALUE.
      VALUE=FX
      CALL MIN(N,1,2,D(1),S,VALUE,.FALSE.,F,X,T,MACHEP,H)
      IF (S.GT.0.D0) GO TO 50
           DO 45 I=1,N
45              V(I,1)=-V(I,1)
50    IF (SF.GT.0.9D0*D(1).AND.0.9D0*SF.LT.D(1)) GO TO 70
           DO 60 I=2,N
60              D(I)=0.D0
C
C.....THE INNER LOOP STARTS HERE.....
70    DO 170 K=2,N
           DO 75 I=1,N
75              Y(I)=X(I)
           SF=FX
           IF (KT.GT.0) ILLC=.TRUE.
80         KL=K
           DF=0.D0
C
C.....A RANDOM STEP FOLLOWS (TO AVOID RESOLUTION VALLEYS).
C     PRAXIS ASSUMES THAT RANDOM RETURNS A RANDOM NUMBER UNIFORMLY
C     DISTRIBUTED IN (0,1).
C
           IF(.NOT.ILLC) GO TO 95
                DO 90 I=1,N
                     S=(0.1D0*LDT+T2*(10**KT))*(RANF(N)-0.5D0)
                     Z(I)=S
                     DO 85 J=1,N
85                        X(J)=X(J)+S*V(J,I)
90              CONTINUE
                FX=F(X,N)
                NF=NF+1
C
C.....MINIMIZE ALONG THE "NON-CONJUGATE" DIRECTIONS V(*,K),...,V(*,N)
C
95         DO 105 K2=K,N
                SL=FX
                S=0.D0
                VALUE=FX
                CALL MIN(N,K2,2,D(K2),S,VALUE,.FALSE.,F,X,T,MACHEP,H)
                IF (ILLC) GO TO 97
                     S=SL-FX
                     GO TO 99
97              S=D(K2)*((S+Z(K2))**2)
99              IF (DF.GT.S) GO TO 105
                     DF=S
                     KL=K2
105        CONTINUE
           IF (ILLC.OR.(DF.GE.DABS((100*MACHEP)*FX))) GO TO 110
C
C.....IF THERE WAS NOT MUCH IMPROVEMENT ON THE FIRST TRY, SET
C     ILLC=TRUE AND START THE INNER LOOP AGAIN.....
C
           ILLC=.TRUE.
           GO TO 80
110        IF (K.EQ.2.AND.PRIN.GT.1) CALL VCPRNT(1,D,N)
C
C.....MINIMIZE ALONG THE "CONJUGATE" DIRECTIONS V(*,1),...,V(*,K-1)
C
           KM1=K-1
           DO 120 K2=1,KM1
           S=0
           VALUE=FX
           CALL MIN(N,K2,2,D(K2),S,VALUE,.FALSE.,F,X,T,MACHEP,H)
120        CONTINUE
           F1=FX
           FX=SF
           LDS=0
           DO 130 I=1,N
                SL=X(I)
                X(I)=Y(I)
                SL=SL-Y(I)
                Y(I)=SL
130             LDS=LDS+SL*SL
           LDS=DSQRT(LDS)
           IF (LDS.LE.SMALL) GO TO 160
C
C.....DISCARD DIRECTION V(*,KL).
C     IF NO RANDOM STEP WAS TAKEN, V(*,KL) IS THE "NON-CONJUGATE"
C     DIRECTION ALONG WHICH THE GREATEST IMPROVEMENT WAS MADE.....
C
           KLMK=KL-K
           IF (KLMK.LT.1) GO TO 141
           DO 140 II=1,KLMK
                I=KL-II
                DO 135 J=1,N
135                  V(J,I+1)=V(J,I)
140             D(I+1)=D(I)
141        D(K)=0
           DO 145 I=1,N
145             V(I,K)=Y(I)/LDS
C
C.....MINIMIZE ALONG THE NEW "CONJUGATE" DIRECTION V(*,K), WHICH IS
C     THE NORMALIZED VECTOR:  (NEW X) - (0LD X).....
C
           VALUE=F1
           CALL MIN(N,K,4,D(K),LDS,VALUE,.TRUE.,F,X,T,MACHEP,H)
           IF (LDS.GT.0.D0) GO TO 160
                LDS=-LDS
                DO 150 I=1,N
150                  V(I,K)=-V(I,K)
160        LDT=LDFAC*LDT
           IF (LDT.LT.LDS) LDT=LDS
           IF (PRIN.GT.0) CALL PRINT(N,X,PRIN,FMIN)
           T2=0.D0
           DO 165 I=1,N
165             T2=T2+X(I)**2
           T2=M2*DSQRT(T2)+T
C
C.....SEE WHETHER THE LENGTH OF THE STEP TAKEN SINCE STARTING THE
C     INNER LOOP EXCEEDS HALF THE TOLERANCE.....
C
           IF (LDT.GT.(0.5*T2)) KT=-1
           KT=KT+1
           IF (KT.GT.KTM) GO TO 400
170   CONTINUE
C.....THE INNER LOOP ENDS HERE.
C
C     TRY QUADRATIC EXTRAPOLATION IN CASE WE ARE IN A CURVED VALLEY.
C
171   CALL QUAD(N,F,X,T,MACHEP,H)
      DN=0.D0
      DO 175 I=1,N
           D(I)=1.D0/DSQRT(D(I))
           IF (DN.LT.D(I)) DN=D(I)
175   CONTINUE
      IF (PRIN.GT.3) CALL MAPRNT(1,V,IDIM,N)
      DO 180 J=1,N
           S=D(J)/DN
           DO 180 I=1,N
180             V(I,J)=S*V(I,J)
C
C.....SCALE THE AXES TO TRY TO REDUCE THE CONDITION NUMBER.....
C
      IF (SCBD.LE.1.D0) GO TO 200
           S=VLARGE
           DO 185 I=1,N
                SL=0.D0
                DO 182 J=1,N
182                  SL=SL+V(I,J)*V(I,J)
                Z(I)=DSQRT(SL)
                IF (Z(I).LT.M4) Z(I)=M4
                IF (S.GT.Z(I)) S=Z(I)
185        CONTINUE
           DO 195 I=1,N
                SL=S/Z(I)
                Z(I)=1.D0/SL
                IF (Z(I).LE.SCBD) GO TO 189
                     SL=1.D0/SCBD
                     Z(I)=SCBD
189             DO 190 J=1,N
190                  V(I,J)=SL*V(I,J)
195        CONTINUE
C
C.....CALCULATE A NEW SET OF ORTHOGONAL DIRECTIONS BEFORE REPEATING
C     THE MAIN LOOP.
C     FIRST TRANSPOSE V FOR MINFIT:
C
200   DO 220 I=2,N
           IM1=I-1
           DO 210 J=1,IM1
                S=V(I,J)
                V(I,J)=V(J,I)
210             V(J,I)=S
220   CONTINUE
C
C.....CALL MINFIT TO FIND THE SINGULAR VALUE DECOMPOSITION OF V.
C     THIS GIVES THE PRINCIPAL VALUES AND PRINCIPAL DIRECTIONS OF THE
C     APPROXIMATING QUADRATIC FORM WITHOUT SQUARING THE CONDITION
C     NUMBER.....
C
      CALL MINFIT(IDIM,N,MACHEP,VSMALL,V,D)
C
C.....UNSCALE THE AXES.....
C
      IF (SCBD.LE.1.D0) GO TO 250
           DO 230 I=1,N
                S=Z(I)
                DO 225 J=1,N
225                  V(I,J)=S*V(I,J)
230        CONTINUE
           DO 245 I=1,N
                S=0.D0
                DO 235 J=1,N
235                  S=S+V(J,I)**2
                S=DSQRT(S)
                D(I)=S*D(I)
                S=1/S
                DO 240 J=1,N
240                  V(J,I)=S*V(J,I)
245        CONTINUE
C
250   DO 270 I=1,N
           DNI=DN*D(I)
           IF (DNI.GT.LARGE) GO TO 265
                IF (DNI.LT.SMALL) GO TO 260
                     D(I)=1/(DNI*DNI)
                     GO TO 270
260             D(I)=VLARGE
                GO TO 270
265        D(I)=VSMALL
270   CONTINUE
C
C.....SORT THE EIGENVALUES AND EIGENVECTORS.....
C
      CALL SORT(IDIM,N,D,V)
      DMIN=D(N)
      IF (DMIN.LT.SMALL) DMIN=SMALL
      ILLC=.FALSE.
      IF (M2*D(1).GT.DMIN) ILLC=.TRUE.
      IF (PRIN.GT.1.AND.SCBD.GT.1.D0) CALL VCPRNT(2,Z,N)
      IF (PRIN.GT.1) CALL VCPRNT(3,D,N)
      IF (PRIN.GT.3) CALL MAPRNT(2,V,IDIM,N)
C.....THE MAIN LOOP ENDS HERE.....
C
      GO TO 40
C
C.....RETURN.....
C
400   IF (PRIN.GT.0) CALL VCPRNT(4,X,N)
      PRAXIS=FX
      RETURN
      END
      SUBROUTINE MINFIT(M,N,MACHEP,TOL,AB,Q)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 MACHEP
      DIMENSION AB(M,N),Q(N)
C...AN IMPROVED VERSION OF MINFIT (SEE GOLUB AND REINSCH, 1969)
C   RESTRICTED TO M=N,P=0.
C   THE SINGULAR VALUES OF THE ARRAY AB ARE RETURNED IN Q AND AB IS
C   OVERWRITTEN WITH THE ORTHOGONAL MATRIX V SUCH THAT U.DIAG(Q) = AB.V,
C   WHERE U IS ANOTHER ORTHOGONAL MATRIX.
      DIMENSION E(100)
C...HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM...
      IF (N.EQ.1) GO TO 200
      EPS = MACHEP
      G = 0.D0
      X = 0.D0
      DO 11 I=1,N
         E(I) = G
         S = 0.D0
         L = I + 1
         DO 1 J=I,N
1           S = S + AB(J,I)**2
         G = 0.D0
         IF (S.LT.TOL) GO TO 4
            F = AB(I,I)
            G = DSQRT(S)
            IF (F.GE.0.D0) G = -G
            H = F*G - S
            AB(I,I)=F-G
            IF (L.GT.N) GO TO 4
            DO 3 J=L,N
               F = 0.D0
               DO 2 K=I,N
2                 F = F + AB(K,I)*AB(K,J)
               F = F/H
               DO 3 K=I,N
3                 AB(K,J) = AB(K,J) + F*AB(K,I)
4        Q(I) = G
         S = 0.D0
         IF (I.EQ.N) GO TO 6
         DO 5 J=L,N
5           S = S + AB(I,J)*AB(I,J)
6        G = 0.D0
         IF (S.LT.TOL) GO TO 10
            IF (I.EQ.N) GO TO 16
            F = AB(I,I+1)
16          G = DSQRT(S)
            IF (F.GE.0.D0) G = -G
            H = F*G - S
            IF (I.EQ.N) GO TO 10
            AB(I,I+1) = F - G
            DO 7 J=L,N
7              E(J) = AB(I,J)/H
            DO 9 J=L,N
               S = 0.D0
               DO 8 K=L,N
8                 S = S + AB(J,K)*AB(I,K)
               DO 9 K=L,N
9                 AB(J,K) = AB(J,K) + S*E(K)
10       Y = DABS(Q(I)) + DABS(E(I))
11       IF (Y.GT.X) X = Y
C...ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS...
      AB(N,N) = 1.D0
      G = E(N)
      L = N
      DO 25 II=2,N
         I = N - II + 1
         IF (G.EQ.0.D0) GO TO 23
         H = AB(I,I+1)*G
         DO 20 J=L,N
20          AB(J,I) = AB(I,J)/H
         DO 22 J=L,N
            S = 0.D0
            DO 21 K=L,N
21             S = S + AB(I,K)*AB(K,J)
            DO 22 K=L,N
22             AB(K,J) = AB(K,J) + S*AB(K,I)
23       DO 24 J=L,N
            AB(I,J) = 0.D0
24          AB(J,I) = 0.D0
         AB(I,I) = 1.D0
         G = E(I)
25       L = I
C...DIAGONALIZATION OF THE BIDIAGONAL FORM...
100   EPS = EPS*X
      DO 150 KK=1,N
         K = N - KK + 1
         KT = 0
101      KT = KT + 1
         IF (KT.LE.30) GO TO 102
            E(K) = 0.D0
            WRITE (6,1000)
1000        FORMAT (' QR FAILED')
102      DO 103 LL2=1,K
            L2 = K - LL2 + 1
            L = L2
            IF (DABS(E(L)).LE.EPS) GO TO 120
            IF (L.EQ.1) GO TO 103
            IF (DABS(Q(L-1)).LE.EPS) GO TO 110
103         CONTINUE
C...CANCELLATION OF E(L) IF L>1...
110      C = 0.D0
         S = 1.D0
         DO 116 I=L,K
            F = S*E(I)
            E(I) = C*E(I)
            IF (DABS(F).LE.EPS) GO TO 120
            G = Q(I)
C...Q(I) = H = DSQRT(G*G + F*F)...
            IF (DABS(F).LT.DABS(G)) GO TO 113
            IF (F) 112,111,112
111         H = 0.D0
            GO TO 114
112         H = DABS(F)*DSQRT(1 + (G/F)**2)
            GO TO 114
113         H = DABS(G)*DSQRT(1 + (F/G)**2)
114         Q(I) = H
            IF (H.NE.0.D0) GO TO 115
               G = 1.D0
               H = 1.D0
115         C = G/H
116         S = -F/H
C...TEST FOR CONVERGENCE...
120      Z = Q(K)
         IF (L.EQ.K) GO TO 140
C...SHIFT FROM BOTTOM 2*2 MINOR...
         X = Q(L)
         Y = Q(K-1)
         G = E(K-1)
         H = E(K)
         F = ((Y - Z)*(Y + Z) + (G - H)*(G + H))/(2*H*Y)
         G = DSQRT(F*F + 1.0D0)
         TEMP = F - G
         IF (F.GE.0.D0) TEMP = F + G
         F = ((X - Z)*(X + Z) + H*(Y/TEMP - H))/X
C...NEXT QR TRANSFORMATION...
         C = 1.D0
         S = 1.D0
         LP1 = L + 1
         IF (LP1.GT.K) GO TO 133
         DO 132 I=LP1,K
            G = E(I)
            Y = Q(I)
            H = S*G
            G = G*C
            IF (DABS(F).LT.DABS(H)) GO TO 123
            IF (F) 122,121,122
121         Z = 0.D0
            GO TO 124
122         Z = DABS(F)*DSQRT(1 + (H/F)**2)
            GO TO 124
123         Z = DABS(H)*DSQRT(1 + (F/H)**2)
124         E(I-1) = Z
            IF (Z.NE.0.D0) GO TO 125
               F = 1.D0
               Z = 1.D0
125         C = F/Z
            S = H/Z
            F = X*C + G*S
            G = -X*S + G*C
            H = Y*S
            Y = Y*C
            DO 126 J=1,N
               X = AB(J,I-1)
               Z = AB(J,I)
               AB(J,I-1) = X*C + Z*S
126            AB(J,I) = -X*S + Z*C
            IF (DABS(F).LT.DABS(H)) GO TO 129
            IF (F) 128,127,128
127         Z = 0.D0
            GO TO 130
128         Z = DABS(F)*DSQRT(1 + (H/F)**2)
            GO TO 130
129         Z = DABS(H)*DSQRT(1 + (F/H)**2)
130         Q(I-1) = Z
            IF (Z.NE.0.D0) GO TO 131
               F = 1.D0
               Z = 1.D0
131         C = F/Z
            S = H/Z
            F = C*G + S*Y
132         X = -S*G + C*Y
133      E(L) = 0.D0
         E(K) = F
         Q(K) = X
         GO TO 101
C...CONVERGENCE:  Q(K) IS MADE NON-NEGATIVE...
140      IF (Z.GE.0.D0) GO TO 150
         Q(K) = -Z
         DO 141 J=1,N
141         AB(J,K) = -AB(J,K)
150      CONTINUE
      RETURN
200   Q(1) = AB(1,1)
      AB(1,1) = 1.D0
      RETURN
      END
      SUBROUTINE MIN(N,J,NITS,D2,X1,F1,FK,F,X,T,MACHEP,H)
      IMPLICIT REAL*8 (A-H,O-Z)
      EXTERNAL F
      LOGICAL FK
      REAL*8 MACHEP,X(N),LDT
      DIMENSION V(100,100),Q0(100),Q1(100)
      COMMON /GLOBAL/ FX,LDT,DMIN,NF,NL
     .       /Q/ V,Q0,Q1,QA,QB,QC,QD0,QD1,QF1
C...THE SUBROUTINE MIN MINIMIZES F FROM X IN THE DIRECTION V(*,J) UNLESS
C   J IS LESS THAN 1, WHEN A QUADRATIC SEARCH IS MADE IN THE PLANE
C   DEFINED BY Q0,Q1,X.
C   D2 IS EITHER ZERO OR AN APPROXIMATION TO HALF F".
C   ON ENTRY, X1 IS AN ESTIMATE OF THE DISTANCE FROM X TO THE MINIMUM
C   ALONG V(*,J) (OR, IF J=0, A CURVE).  ON RETURN, X1 IS THE DISTANCE
C   FOUND.
C   IF FK=.TRUE., THEN F1 IS FLIN(X1).  OTHERWISE X1 AND F1 ARE IGNORED
C   ON ENTRY UNLESS FINAL FX IS GREATER THAN F1.
C   NITS CONTROLS THE NUMBER OF TIMES AN ATTEMPT WILL BE MADE TO HALVE
C   THE INTERVAL.
      LOGICAL DZ
      REAL*8 M2,M4
      SMALL = MACHEP**2
      M2 = DSQRT(MACHEP)
      M4 = DSQRT(M2)
      SF1 = F1
      SX1 = X1
      K = 0
      XM = 0.D0
      FM = FX
      F0 = FX
      DZ = D2.LT.MACHEP
C...FIND THE STEP SIZE...
      S = 0.D0
      DO 1 I=1,N
1        S = S + X(I)**2
      S = DSQRT(S)
      TEMP = D2
      IF (DZ) TEMP = DMIN
      T2 = M4*DSQRT(DABS(FX)/TEMP + S*LDT) + M2*LDT
      S = M4*S + T
      IF (DZ.AND.T2.GT.S) T2 = S
      T2 = DMAX1(T2,SMALL)
      T2 = DMIN1(T2,.01D0*H)
      IF (.NOT.FK.OR.F1.GT.FM) GO TO 2
      XM = X1
      FM = F1
2     IF (FK.AND.DABS(X1).GE.T2) GO TO 3
      TEMP=1.D0
      IF (X1.LT.0.D0) TEMP=-1.D0
      X1=TEMP*T2
      F1 = FLIN(N,J,X1,F,X,NF)
3     IF (F1.GT.FM) GO TO 4
      XM = X1
      FM = F1
4     IF (.NOT.DZ) GO TO 6
C...EVALUATE FLIN AT ANOTHER POINT AND ESTIMATE THE SECOND DERIVATIVE...
      X2 = -X1
      IF (F0.GE.F1) X2 = 2.D0*X1
      F2 = FLIN(N,J,X2,F,X,NF)
      IF (F2.GT.FM) GO TO 5
         XM = X2
         FM = F2
5     D2 = (X2*(F1 - F0)-X1*(F2 - F0))/((X1*X2)*(X1 - X2))
C...ESTIMATE THE FIRST DERIVATIVE AT 0...
6     D1 = (F1 - F0)/X1 - X1*D2
      DZ = .TRUE.
C...PREDICT THE MINIMUM...
      IF (D2.GT.SMALL) GO TO 7
         X2 = H
         IF (D1.GE.0.D0) X2 = -X2
         GO TO 8
7        X2 = (-.5D0*D1)/D2
8     IF (DABS(X2).LE.H) GO TO 11
         IF (X2) 9,9,10
9        X2 = -H
         GO TO 11
10       X2 = H
C...EVALUATE F AT THE PREDICTED MINIMUM...
11    F2 = FLIN(N,J,X2,F,X,NF)
      IF (K.GE.NITS.OR.F2.LE.F0) GO TO 12
C...NO SUCCESS, SO TRY AGAIN...
         K = K + 1
         IF (F0.LT.F1.AND.(X1*X2).GT.0.D0) GO TO 4
         X2 = 0.5D0*X2
         GO TO 11
C...INCREMENT THE ONE-DIMENSIONAL SEARCH COUNTER...
12    NL = NL + 1
      IF (F2.LE.FM) GO TO 13
      X2 = XM
      GO TO 14
13    FM = F2
C...GET A NEW ESTIMATE OF THE SECOND DERIVATIVE...
14    IF (DABS(X2*(X2 - X1)).LE.SMALL) GO TO 15
         D2 = (X2*(F1-F0) - X1*(FM-F0))/((X1*X2)*(X1 - X2))
         GO TO 16
15       IF (K.GT.0) D2 = 0.D0
16    IF (D2.LE.SMALL) D2 = SMALL
      X1 = X2
      FX = FM
      IF (SF1.GE.FX) GO TO 17
         FX = SF1
         X1 = SX1
C...UPDATE X FOR LINEAR BUT NOT PARABOLIC SEARCH...
17    IF (J.EQ.0) RETURN
      DO 18 I=1,N
18       X(I) = X(I) + X1*V(I,J)
      RETURN
      END
      REAL*8 FUNCTION FLIN (N,J,L,F,X,NF)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 L,X(N)
      DIMENSION V(100,100),Q0(100),Q1(100)
C...FLIN IS THE FUNCTION OF ONE REAL VARIABLE L THAT IS MINIMIZED
C   BY THE SUBROUTINE MIN...
      COMMON /Q/ V,Q0,Q1,QA,QB,QC,QD0,QD1,QF1
      DIMENSION T(100)
      IF (J .EQ. 0) GO TO 2
C...THE SEARCH IS LINEAR...
      DO 1 I=1,N
1        T(I) = X(I) + L*V(I,J)
      GO TO 4
C...THE SEARCH IS ALONG A PARABOLIC SPACE CURVE...
2     QA = (L*(L - QD1))/(QD0*(QD0 + QD1))
      QB = ((L + QD0)*(QD1 - L))/(QD0*QD1)
      QC = (L*(L + QD0))/(QD1*(QD0 + QD1))
      DO 3 I=1,N
3        T(I) = (QA*Q0(I) + QB*X(I)) + QC*Q1(I)
C...THE FUNCTION EVALUATION COUNTER NF IS INCREMENTED...
4     NF = NF + 1
      FLIN = F(T,N)
      RETURN
      END
      SUBROUTINE SORT(M,N,D,V)
      REAL*8 D(N),V(M,N)
C...SORTS THE ELEMENTS OF D(N) INTO DESCENDING ORDER AND MOVES THE
C   CORRESPONDING COLUMNS OF V(N,N).
C   M IS THE ROW DIMENSION OF V AS DECLARED IN THE CALLING PROGRAM.
      REAL*8 S
      IF (N.EQ.1) RETURN
      NM1 = N - 1
      DO 3 I = 1,NM1
         K=I
         S = D(I)
         IP1 = I + 1
         DO 1 J = IP1,N
            IF (D(J) .LE. S) GO TO 1
            K = J
            S = D(J)
1           CONTINUE
         IF (K .LE. I) GO TO 3
         D(K) = D(I)
         D(I) = S
         DO 2 J = 1,N
            S = V(J,I)
            V(J,I) = V(J,K)
2           V(J,K) = S
3        CONTINUE
      RETURN
      END
      SUBROUTINE QUAD(N,F,X,T,MACHEP,H)
      IMPLICIT REAL*8 (A-H,O-Z)
      EXTERNAL F
C...QUAD LOOKS FOR THE MINIMUM OF F ALONG A CURVE DEFINED BY Q0,Q1,X...
      REAL*8 X(N),MACHEP,LDT,L
      DIMENSION V(100,100),Q0(100),Q1(100)
      COMMON /GLOBAL/ FX,LDT,DMIN,NF,NL
     .       /Q/ V,Q0,Q1,QA,QB,QC,QD0,QD1,QF1
      S = FX
      FX = QF1
      QF1 = S
      QD1 = 0.D0
      DO 1 I=1,N
         S = X(I)
         L = Q1(I)
         X(I) = L
         Q1(I) = S
1        QD1 = QD1 + (S-L)**2
      QD1 = DSQRT(QD1)
      L = QD1
      S = 0.D0
      IF (QD0. LE. 0.D0 .OR. QD1 .LE. 0.D0 .OR. NL .LT. 3*N*N) GO TO 2
      VALUE=QF1
      CALL MIN(N,0,2,S,L,VALUE,.TRUE.,F,X,T,MACHEP,H)
      QA = (L*(L-QD1))/(QD0*(QD0+QD1))
      QB = ((L+QD0)*(QD1-L))/(QD0*QD1)
      QC = (L*(L+QD0))/(QD1*(QD0+QD1))
      GO TO 3
2     FX = QF1
      QA = 0.D0
      QB = QA
      QC = 1.D0
3     QD0 = QD1
      DO 4 I=1,N
         S = Q0(I)
         Q0(I) = X(I)
4        X(I) = (QA*S + QB*X(I)) + QC*Q1(I)
      RETURN
      END
      SUBROUTINE VCPRNT(OPTION,V,N)
      REAL*8 V(N)
      INTEGER OPTION
      GO TO (1,2,3,4),OPTION
1     WRITE (6,101) (V(I),I=1,N)
      RETURN
2     WRITE (6,102) (V(I),I=1,N)
      RETURN
3     WRITE (6,103) (V(I),I=1,N)
      RETURN
4     WRITE (6,104) (V(I),I=1,N)
      RETURN
101   FORMAT (/' THE SECOND DIFFERENCE ARRAY D(*) IS:'/
     .        (E32.14,4E25.14))
102   FORMAT (/' THE SCALE FACTORS ARE:'/(E32.14,4E25.14))
103   FORMAT (/' THE APPROXIMATING QUADRATIC FORM HAS THE PRINCIPAL VALU
     .ES:'/(E32.14,4E25.14))
104   FORMAT (/' X IS:',E26.14/(E32.14))
      END
      SUBROUTINE PRINT(N,X,PRIN,FMIN)
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER PRIN
      REAL*8 X(N),LN,LDT
      COMMON /GLOBAL/ FX,LDT,DMIN,NF,NL
      WRITE (6,101) NL,NF,FX
      IF (FX.LE.FMIN) GO TO 1
      LN = DLOG10(FX-FMIN)
      WRITE (6,102) FMIN,LN
      GO TO 2
1     WRITE (6,103) FMIN
2     IF (N.GT.4.AND.PRIN.LE.2) RETURN
      WRITE (6,104) (X(I),I=1,N)
      RETURN
101   FORMAT (/' AFTER',I6,
     . ' LINEAR SEARCHES, THE FUNCTION HAS BEEN EVALUATED',I6,
     . ' TIMES.  THE SMALLEST VALUE FOUND IS F(X) = ',E21.14)
102   FORMAT (' LOG (F(X)-',E21.14,') = ',E21.14)
103   FORMAT (' LOG (F(X)-',E21.14,') IS UNDEFINED.')
104   FORMAT (' X IS:',E26.14/(E32.14))
      END
      SUBROUTINE MAPRNT(OPTION,V,M,N)
      REAL*8 V(M,N)
      INTEGER OPTION,LOW,UPP
C...THE SUBROUTINE MAPRNT PRINTS THE COLUMNS OF THE NXN MATRIX V
C   WITH A HEADING AS SPECIFIED BY OPTION.
C   M IS THE ROW DIMENSION OF V AS DECLARED IN THE CALLING PROGRAM...
      LOW = 1
      UPP = 5
      GO TO (1,2),OPTION
1     WRITE (6,101)
101   FORMAT (/' THE NEW DIRECTIONS ARE:')
      GO TO 3
2     WRITE (6,102)
102   FORMAT (' AND THE PRINCIPAL AXES:')
3     IF (N.LT.UPP) UPP = N
      DO 4 I=1,N
4        WRITE (6,104) (V(I,J),J=LOW,UPP)
      LOW = LOW + 5
      IF (N.LT.LOW) RETURN
      UPP = UPP + 5
      WRITE (6,103)
      GO TO 3
103   FORMAT (' ')
104   FORMAT (E32.14,4E25.14)
      END
      REAL*8 FUNCTION RANDOM(NAUGHT)
      REAL*8 RAN1,RAN3(127),HALF
      INTEGER RAN2,Q,R
      LOGICAL INIT
      DATA INIT/.FALSE./
      IF (INIT) GO TO 3
      R = MOD(NAUGHT,8190) + 1
      RAN2 = 128
      DO 2 I=1,127
         RAN2 = RAN2 - 1
         RAN1 = -2.D0**55
         DO 1 J=1,7
            R = MOD(1756*R,8191)
            Q = R/32
1           RAN1 = (RAN1 + Q)*(1.0D0/256)
2        RAN3(RAN2) = RAN1
      INIT = .TRUE.
3     IF (RAN2.EQ.1) RAN2 = 128
      RAN2 = RAN2 - 1
      RAN1 = RAN1 + RAN3(RAN2)
      HALF = .5D0
      IF (RAN1.GE.0.D0) HALF = -HALF
      RAN1 = RAN1 + HALF
      RAN3(RAN2) = RAN1
      RANDOM = RAN1 + .5D0
      RETURN
      END

      REAL*8 FUNCTION RANF(NAUGHT)
      DATA IA/16807/,IC/2147483647/,IQ/127773/,IR/2836/
      COMMON /CSEED/ ISEED
        IH = ISEED/IQ
        IL = MOD(ISEED,IQ)
        IT = IA*IL-IR*IH
        IF(IT.GT.0) THEN
          ISEED = IT
        ELSE
          ISEED = IC+IT
        END IF
        RANF = ISEED/FLOAT(IC)
      RETURN
      END

